##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(ArchR)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){

    ## 所有的细胞的,计算maxdelt
    comine_data_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/qc_atac_v3/germ/"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################
## 导入数据
b <- load(comine_data_file)
## testis_combined_peak_combineRNA

###########################################################################################
## ATAC的开放矩阵
GSM_se <- getMatrixFromProject(testis_combined_peak_combineRNA, useMatrix="GeneScoreMatrix")
GSM_mat <- assays(GSM_se)$GeneScoreMatrix
rownames(GSM_mat) <- rowData(GSM_se)$name

## 转化为数值矩阵
GSM_mat_num <- apply(GSM_mat , 1 , as.numeric)
rownames(GSM_mat_num) <- colnames(GSM_mat)
GSM_mat_num <- t(GSM_mat_num)
GSM_mat_num <- GSM_mat_num[,rownames(testis_combined_peak_combineRNA@cellColData)]

## 构造atac矩阵
sco_atac <- sapply(unique(unique(testis_combined_peak_combineRNA@cellColData$cell_type)),function(x){
    print(x)
    sapply(unique(rownames(GSM_mat_num)),function(y){
        mean(
            as.numeric(
                as.vector(GSM_mat_num[y,which(testis_combined_peak_combineRNA@cellColData$cell_type==x)])
            )
        )
    })
})


sco_exp2 <- data.frame(sco_atac)
sco_exp2$gene <- rownames(sco_exp2)
sco_exp2 <- sco_exp2[,c( ncol(sco_exp2) , (1:ncol(sco_exp2)-1) )]

out_file <- paste0( out_path , "/GeneScore.MeanByCellType.tsv" )
write.table( sco_exp2 , out_file , row.names = F , sep = "\t" )